In [1]:
#
import numpy as np
import os, sys

sys.path.insert(0, '/global/homes/q/qmxp55/DESI/bgstargets/py')

from io_ import get_sweep_whole, getBGSbits, flux_to_mag
from io_ import get_random, get_isdesi, get_dict, bgsmask, get_reg, get_svfields, get_svfields_fg, gaiaAEN
from cuts import getGeoCuts, bgsbut
from QA import getStats, flow, mollweide, mycmap, plot_sysdens, overdensity, hexbin
from postages_images import postages_circle

import healpy as hp
import astropy.io.fits as fits
import fitsio
import matplotlib.pyplot as plt
from matplotlib import gridspec
import matplotlib
from astropy.coordinates import SkyCoord
import astropy.units as units

%load_ext autoreload
%autoreload 2

import warnings
warnings.filterwarnings('ignore')
In [2]:
#
dr = 'dr8'
survey = 'north' #is either south (DECaLS+DES) or north (BASS/MzLS)
version = '0.1.0'
filesdir = '/global/cscratch1/sd/qmxp55/bgstargets_output/'
Nranfiles = 3
reg = 'svfields_fg'
dec_resol_ns = 32.375
if (dr[:3] == 'dr9'): Nranfiles = 20 # because the randoms for dr9d have a density of 100000 = 5000*20

pathdir = os.path.abspath(os.getcwd())+'/%s-%s_%s_%s_draft_plots' %(dr, survey, reg, version)
ispathdir = os.path.isdir(pathdir)
if not ispathdir: os.mkdir(pathdir)
    
if reg[:8] == 'svfields': desifootprint = False
else: desifootprint = True

# for healpy
dr8pix     = '/project/projectdirs/desi/target/catalogs/dr8/0.31.1/pixweight/pixweight-dr8-0.31.1.fits'
hdr          = fits.getheader(dr8pix,1)
nside,nest   = hdr['hpxnside'],hdr['hpxnest']
npix         = hp.nside2npix(nside)
pixarea      = hp.nside2pixarea(nside,degrees=True)
In [3]:
#load catalogues
if dr[:3] == 'dr9': N = 1
else: N = Nranfiles
    
cat = np.load(filesdir+dr+'/'+version+'/'+'bgstargets-'+survey+'.npy')
cat_ex = np.load(filesdir+dr+'/'+version+'/'+'extra-'+survey+'_n256.npy')
ran = np.load(filesdir+dr+'/'+dr+'_random_N'+str(N)+'.npy')
ran_ex = np.load(filesdir+dr+'/'+'extra_random_N'+str(N)+'_n256.npy')
In [17]:
hppix_ran = ran_ex['hppix']
ranindesi = ran_ex['desi']
#if svfields don't use desifootprint
#tunning in south as dr9f/dr9g does not entirely cover southern ages region
if (reg == 'svfields_fg') & (survey == 'south'): tun, tun_ran = (cat['DEC'] < 30), (ran['DEC'] < 30)
else: tun, tun_ran = np.ones_like(cat['RA'], dtype=bool), np.ones_like(ran['RA'], dtype=bool)
#tunning in north as dr9f/dr9g does not entirely cover highstardens_n region
if (reg == 'svfields_fg') & (survey == 'north'): tun2, tun2_ran = (cat['RA'] < 280), (ran['RA'] < 280)
else: tun2, tun2_ran = np.ones_like(cat['RA'], dtype=bool), np.ones_like(ran['RA'], dtype=bool)
    
if not desifootprint: raninreg = (ran_ex[reg]) & (ran_ex[survey]) & (tun_ran) & (tun2_ran)
else: raninreg = (ran_ex[reg]) & (ran_ex['desi']) & (ran_ex[survey]) & (tun_ran) & (tun2_ran)

hppix_cat = cat_ex['hppix']
catindesi = cat_ex['desi']
#if svfields don't use desifootprintreal amateurreal amateur
if not desifootprint: catinreg = (cat_ex[reg]) & (cat_ex[survey]) & (tun) & (tun2)
else: catinreg = catinreg = (cat_ex[reg]) & (cat_ex['desi']) & (cat_ex[survey]) & (tun) & (tun2)

Filter catalogue and randoms by region

In [18]:
#RANDOMS
keep = raninreg
ran = ran[keep]
ran_ex = ran_ex[keep]
ranindesi = ranindesi[keep]
hppix_ran = hppix_ran[keep]
raninreg = raninreg[keep]

#CATALOGUE
keepC = catinreg
cat = cat[keepC]
cat_ex = cat_ex[keepC]
catindesi = catindesi[keepC]
hppix_cat = hppix_cat[keepC]
catinreg = catinreg[keepC]
In [19]:
#bgsmask = bgsmask()
rancuts = getGeoCuts(ran, randoms=True)

getting area of reg with healpix pixels (for flow charts...)

NOTE: if using SVFIELDS footprint, use desifootprint=False otherwise the area will be less than the current SVFIELDS due desi footprint chops out part of it.

In [21]:
#Notes: This hpdict is used to get the area only using the randoms. 
#This use the randoms within the DESI footprint and without any masking (maskrand=None) as we want the area without wholes.
hpdict0 = get_dict(cat=None, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=None, 
                      maskrand=None, maskcat=None, Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None, 
                            desifootprint=desifootprint, namesels=None, target_outputs=False, log=True)

positions and desifotprint DONE...
systematics DONE...
randdens =  15000  ; len randoms =  1665831
bgsfracarea DONE...
regions DONE...
area_all = 111 deg2
bgsarea_south = 0 deg2
bgsarea_decals = 0 deg2
bgsarea_des = 0 deg2
bgsarea_north = 111 deg2
bgsarea_south_n = 0 deg2
bgsarea_south_s = 0 deg2
bgsarea_svfields = 110 deg2
bgsarea_svfields_n = 110 deg2
bgsarea_svfields_s = 0 deg2
bgsarea_svfields_fg = 110 deg2
bgsarea_svfields_fg_n = 110 deg2
bgsarea_svfields_fg_s = 0 deg2
areas DONE...

Healpy pixels catalogue with BGS target density (for dens. sky maps and systematics plots...)

In [22]:
#from astropy.coordinates import SkyCoord
#import astropy.units as units

#c = SkyCoord(cat['RA']*units.degree,cat['DEC']*units.degree, frame='icrs')
#galb = c.galactic.b.value # galb coordinate
In [23]:
#dic with default BGS selection and in DESI footprint
#Appy or not apply LG mask in randoms (rancuts['LG'])?
hpdict = get_dict(cat=cat, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=hppix_cat, 
                      maskrand=((rancuts['BS']) & (rancuts['GC']) & (rancuts['nobs'])),
                          maskcat=None, 
                             Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None, 
                                desifootprint=desifootprint, namesels={'bgs_any':20, 'bright':21, 'faint':22}, 
                                  galb=cat_ex['b'], log=False, survey='bgs')

BGS outline

In [24]:
#
bgs_any = ((cat['BGSBITS'] & 2**(20)) != 0)
bgs_b = ((cat['BGSBITS'] & 2**(21)) != 0)
bgs_f = ((cat['BGSBITS'] & 2**(22)) != 0)

print('mean density of BGS: \t %.2f' %(hpdict['meandens_bgs_any_'+survey]))
#print('total density of BGS: \t %.2f' %(hpdict['dens_bgs_any_'+survey]))
print('Area within BGS: \t %.2f' %(hpdict['bgsarea_'+survey]))

print('No BGS galaxies in Bright and Faint: \t %i \t %i \t %i' %(
    np.sum(bgs_b),
    np.sum(bgs_f),
    np.sum(bgs_any)
))


print('Overall density of BGS: \t %d \t %d \t %d' %(
    np.round(hpdict['dens_bright_'+survey], 0),
    np.round(hpdict['dens_faint_'+survey], 0),
    np.round(hpdict['dens_bgs_any_'+survey], 0)
))

print('morph \t BRIGHT \t FAINT \t OVERALL')
print('------------------------------------')
for morph in list(set(cat['TYPE'][:10000])):
    
    keep = cat['TYPE'] == morph
    print('%s \t %d \t %d \t %d' %(
        morph, 
        np.round(np.sum((keep) & (bgs_b))/hpdict['bgsarea_'+survey], 0),
        np.round(np.sum((keep) & (bgs_f))/hpdict['bgsarea_'+survey], 0),
        np.round(np.sum((keep) & (bgs_any))/hpdict['bgsarea_'+survey], 0),
        ))

mean density of BGS: 	 1393.46
Area within BGS: 	 107.64
No BGS galaxies in Bright and Faint: 	 88205 	 62577 	 150782
Overall density of BGS: 	 819 	 581 	 1401
morph 	 BRIGHT 	 FAINT 	 OVERALL
------------------------------------
REX  	 103 	 146 	 249
EXP  	 250 	 206 	 456
PSF  	 3 	 2 	 5
COMP 	 23 	 3 	 26
DEV  	 440 	 224 	 665
In [25]:
if dr[:3] == 'dr9': 
    morphos = ['REX', 'EXP', 'DEV', 'SER']
    PSF = ['PSF']
else: 
    morphos = ['REX ', 'EXP ', 'DEV ', 'COMP']
    PSF = ['PSF ']
In [26]:
set(cat['TYPE'][:10000])
Out[26]:
{'COMP', 'DEV ', 'EXP ', 'PSF ', 'REX '}

BGS dens. sky maps

In [27]:
from QA import mollweide, mycmap, plot_sysdens
import matplotlib.pyplot as plt
from matplotlib import gridspec
import matplotlib
In [28]:
#
if reg[:8] == 'svfields': reg_ = reg+'_'+survey[0]
else: reg_ = reg
    
fig = plt.figure(figsize=(15, 15))
gs = gridspec.GridSpec(2,1)
    
org          = 120  # centre ra for mollweide plots
projection   = 'mollweide'
cm = mycmap(matplotlib.cm.jet, 10,0,1)
cmr= mycmap(matplotlib.cm.jet_r,10,0,1)

mollweide(hpdict=hpdict, namesel='bright', reg=reg_, projection=projection, n=0, org=org, cm=cm, fig=fig, gs=gs, cval=None, desifootprint=False)
mollweide(hpdict=hpdict, namesel='faint', reg=reg_, projection=projection, n=1, org=org, cm=cm, fig=fig, gs=gs, cval=None, desifootprint=False)
#mollweide(hpdict=hpdict, C=None ,namesel='any', reg='all', projection=projection, n=1, org=org, cm=cm, 
#          fig=fig, ws=ws, perc=(0.3,99.8), title='After linear weights', cval=(84, 2750))

file = '%s/skydens_%s_%s' %(pathdir, dr, reg)
fig.savefig(file+'.png', bbox_inches = 'tight', pad_inches = 0)

Flow charts: Nominal & Galaxy View

In [29]:
#use below code to get target densities of keep and rejected at a particular stage in flowchart...
#use hpdict0 if density computed over the total area (before applyin spatial cuts in randoms)
#use hpdict if density computed over the reduced area (after applying spatial cuts in randoms)
if False:
    t = getStats(cat=cat, hpdict=hpdict0, bgsmask=bgsmask(), rancuts=rancuts, CurrentMask=['QC_FM', 'QC_FI', 'QC_FF'], 
             PrevMask=['BS', 'LG', 'GC', 'nobs', 'SG', 'FMC2', 'CC'], 
                 reg=reg, regcat=catinreg, regran=raninreg)
In [30]:
if reg == 'south': surveylab = 'DECaLS+DES'
elif reg == 'north': surveylab = 'BASS-MzLS'
elif reg == 'desi': surveylab = 'DESI'
else: surveylab = reg
    
flowTitle = dr+'_'+surveylab

#this is important to read the area of the proper region in hpdict0
region = reg
In [31]:
order = [['BS', 'LG', 'GC'], ['nobs'], ['SG'], ['FMC2', 'CC'], ['QC_FM', 'QC_FI', 'QC_FF']]


flowNominal, _, _ = flow(cat=cat, hpdict=hpdict0, bgsmask=bgsmask(), rancuts=rancuts, order=order, reg=region, 
             regcat=catinreg, regran=raninreg, file='%s/flow_main_nominal_%s' %(pathdir, flowTitle), dr=flowTitle, program='main')

flowNominal

Previous Cuts: (None)
Current Cuts: (BS|LG|GC)
Previous Cuts: (BS|LG|GC)
Current Cuts: (nobs)
Previous Cuts: (BS|LG|GC|nobs)
Current Cuts: (SG)
Previous Cuts: (BS|LG|GC|nobs|SG)
Current Cuts: (FMC2|CC)
Previous Cuts: (BS|LG|GC|nobs|SG|FMC2|CC)
Current Cuts: (QC_FM|QC_FI|QC_FF)
Out[31]:
In [32]:
order = [ ['SG'], ['BS', 'LG', 'GC'], ['nobs'], ['QC_FM', 'QC_FI', 'QC_FF'], ['FMC2', 'CC']]

flowGalview, _, _ = flow(cat=cat, hpdict=hpdict0, bgsmask=bgsmask(), rancuts=rancuts, order=order, reg=region, 
             regcat=catinreg, regran=raninreg, file='%s/flow_main_galview_%s' %(pathdir,flowTitle), dr=flowTitle, program='main')

flowGalview

Previous Cuts: (None)
Current Cuts: (SG)
Previous Cuts: (SG)
Current Cuts: (BS|LG|GC)
Previous Cuts: (SG|BS|LG|GC)
Current Cuts: (nobs)
Previous Cuts: (SG|BS|LG|GC|nobs)
Current Cuts: (QC_FM|QC_FI|QC_FF)
Previous Cuts: (SG|BS|LG|GC|nobs|QC_FM|QC_FI|QC_FF)
Current Cuts: (FMC2|CC)
Out[32]:

2D-stacks around Bright Stars and Medium Stars

Firs we do the stacking around the bright stars (BS mask) using BGS but BS.

For the second stacking, we use BGS and do the stacking around the Medium Stars (MS mask).

In [33]:
#star catalogue
stars = np.load('/global/cscratch1/sd/qmxp55/pauline/stars_GAIA_TYCHO_13.npy')
In [34]:
hppix_stars = hp.ang2pix(nside,(90.-stars['DEC'])*np.pi/180.,stars['RA']*np.pi/180.,nest=True) # catalogue hp pixels array
if reg[:8] == 'svfields': starsinreg = get_svfields_fg(stars['RA'], stars['DEC'])
else: starsinreg = get_reg(reg=reg, hppix=hppix_stars)
#cstars = SkyCoord(stars['RA']*units.degree,stars['DEC']*units.degree, frame='icrs')
#galb_stars = cstars.galactic.b.value # galb coordinate
stars = stars[(starsinreg) & (stars['DEC'] > -25)]
In [35]:
#
#bgslist = ['LG', 'GC', 'nobs', 'SG', 'FMC2', 'CC', 'QC_FM', 'QC_FI', 'QC_FF']
    
# get BGS objects without BS bit set on
#bgsbutBS = np.ones_like(cat['RA'], dtype=bool)

#for key in bgslist:
#    keep = ((cat['BGSBITS'] & 2**(bgsmask()[key])) != 0)
#    bgsbutBS &= keep
#bgsbutBS &= cat['RMAG'] < 20
    
bgsbutBS = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['BS'], bgsmask=bgsmask(), rlimit=20)
    
log = True
nbins = 150
    
#Dustin_radii's
mag = np.linspace(0, 20, 50)
BS_radii = []
new_BS_radii = []
for i,j in enumerate(mag):
    BS_radii.append([j, np.minimum(1800., 150. * 2.5**((11. - j)/3.)) * 0.262/1])
    new_BS_radii.append([j, 0.5 * np.minimum(1800., 150. * 2.5**((11. - j)/3.)) * 0.262/1])
    #radius = np.minimum(1800., 150. * 2.5**((11. - mag)/3.)) * 0.262/3600.
        
plt.figure(figsize=(20, 20))
plt.title(r'%s %s' %(dr, survey), size=20)
plt.axis('off')
_ = overdensity(cat[bgsbutBS], stars, BS_radii, 'MAG', 35, density=False, 
                        magbins=[10,13,4], radii_2=new_BS_radii, grid=[2,2], SR=[0., 250.], scaling=False, nbins=nbins, 
                               SR_scaling=3,logDenRat=[-5, 5], radii_bestfit=False, annulus=None, bintype='1', 
                                   filename='%s/2d_stack_BS_%s_%s' %(pathdir, dr, survey), log=log)

mag_bins_len: 4
MAG < 10.00
mag_mean 8.997807343281645
mask_radius2 36.231417202811464
mag_radii MAX: 220.01593742837198 mag_radii MIN: 53.556070423289
mag MAX: 9.993000030517578 mag MIN: 5.363476753234863
density cumu (min, max): (-inf, 86)
density non-cumu (min, max): (-inf, 104)
7251 of inf in density ratio out of a total of 16977
Minimum density ratio = -5, Maximum density ratio = 4.11772
----------------
10.00 < MAG < 11.00
mag_mean 10.56915281087512
mask_radius2 22.430672920527783
mag_radii MAX: 53.38328665799261 mag_radii MIN: 39.32438202738072
mag MAX: 10.999266624450684 mag MIN: 10.003599166870117
density cumu (min, max): (-inf, 36.6)
density non-cumu (min, max): (-76.4, 85.5)
5084 of inf in density ratio out of a total of 16977
Minimum density ratio = -5, Maximum density ratio = 3.3659
----------------
11.00 < MAG < 12.00
mag_mean 11.556748714265565
mask_radius2 16.604552662787317
mag_radii MAX: 39.30752904852582 mag_radii MIN: 29.01030303169343
mag MAX: 11.999972343444824 mag MIN: 11.000593185424805
density cumu (min, max): (-inf, -0.358)
density non-cumu (min, max): (-91, 7.4)
1302 of inf in density ratio out of a total of 16977
Minimum density ratio = -5, Maximum density ratio = 2.29633
----------------
12.00 < MAG < 13.00
mag_mean 12.557205164876324
mask_radius2 12.22993629999663
mag_radii MAX: 29.005276801615953 mag_radii MIN: 21.357443760921427
mag MAX: 12.999948501586914 mag MIN: 12.000547409057617
density cumu (min, max): (-inf, 0.214)
density non-cumu (min, max): (-224, 5.11)
89 of inf in density ratio out of a total of 16977
Minimum density ratio = -5, Maximum density ratio = 1.52721
----------------
In [36]:
plt.figure(figsize=(10,5))

plt.plot(np.transpose(BS_radii)[0], np.transpose(BS_radii)[1], lw=2, c='k', label='BS radii')
plt.plot(np.transpose(new_BS_radii)[0], np.transpose(new_BS_radii)[1], lw=2, c='r', label='NEW BS radii')

plt.legend()
plt.grid()
In [37]:
from QA import circular_mask_radii_func
from io_ import query_catalog_mask
In [38]:
BS_t = query_catalog_mask(cat['RA'], cat['DEC'], stars, BS_radii, nameMag='MAG', diff_spikes=False, 
                             length_radii=None, widht_radii=None, return_diagnostics=True, bestfit=False, log=False)

Total run time: 1.982342 sec
In [39]:
new_BS_t = query_catalog_mask(cat['RA'], cat['DEC'], stars, new_BS_radii, nameMag='MAG', diff_spikes=False, 
                             length_radii=None, widht_radii=None, return_diagnostics=True, bestfit=False, log=False)

Total run time: 2.001950 sec
In [40]:
#masking radius for each object cloosest stars in cat
mask_radii = circular_mask_radii_func(BS_t[1]['w1_source'], BS_radii, bestfit=False)
#masking radius reescaled for each object
mask_radii_res = BS_t[1]['d2d_source']/mask_radii
In [41]:
BS_diff = (BS_t[0]) & (~new_BS_t[0]) & (cat['RMAG'] < 20)
PSF = ((cat['TYPE'] == 'PSF') | (cat['TYPE'] == 'PSF ')) 
refcat = cat['REF_CAT']
if isinstance(np.atleast_1d(refcat)[0], str):
    LG = [(rc[0] == "L") if len(rc) > 0 else False for rc in refcat]
else:
    LG = [(rc.decode()[0] == "L") if len(rc) > 0 else False for rc in refcat]

print('objects within BS radii defs & r<20: \t %i (PSF) \t %i (no-PSF) \t %i (LSLGA)' 
      %(np.sum((BS_diff) & (PSF)), np.sum((BS_diff) & (~PSF)), np.sum((BS_diff) & (LG))))
print('objects within BS radii defs & r<20 & BGS: \t %i (PSF) \t %i (no-PSF) \t %i (LSLGA)' 
      %(np.sum((BS_diff) & (PSF) & (bgsbutBS)), np.sum((BS_diff) & (~PSF) & (bgsbutBS)), np.sum((BS_diff) & (LG))))

objects within BS radii defs & r<20: 	 15716 (PSF) 	 100 (no-PSF) 	 17 (LSLGA)
objects within BS radii defs & r<20 & BGS: 	 1163 (PSF) 	 14 (no-PSF) 	 17 (LSLGA)

$\Delta \vec{r}/r_{BS}$

In [42]:
#
rows, cols = 1, 2
coord = {'r':cat['RMAG'], '$\Delta r/r_{BS}$':mask_radii_res}
area = hpdict0['bgsarea_'+survey]

#debug = cat['RMAG'] < 16
masks = {'BS_old not in BS_new & r<20':(BS_diff), 'BS_old not in BS_new & r<20 & BGS':(BS_diff) & (bgsbutBS)}
hline, vline = 0.5, None
vmin, vmax = 1, None
    
fig    = plt.figure(figsize=(8*cols,6*rows))
gs     = gridspec.GridSpec(rows, cols, hspace=0.25, wspace=0.10)
cmap = plt.get_cmap('coolwarm', 8)
    
for i, key, val in zip(range(len(masks)), masks.keys(), masks.values()):
        
    if (i%cols==0): ylab=True
    else: ylab = False
        
    hexbin(coord=coord, catmask=(val), n=i, bins='log', title=key, cmap=cmap, 
               ylab=ylab, xlab=True, vline=vline, hline=hline, fig=fig, gs=gs, vmin=vmin, vmax=vmax, mincnt=1, fmcline=False, 
                   file='%s/r_mask_radii_BSold_not_BSnew_%s' %(pathdir, survey), fracs=False, area=area, cbar='horizontal',
                      xlim=[12, 20.1], ylim=[0.4, 1.0])
In [43]:
#
veto = {'BGS':(bgs_any),'Not BGS':(~bgs_any)}

info = {'r':cat['RMAG'], 'T':cat['TYPE'].astype(str), 'ref':cat['REF_CAT'].astype(str), 'd2d':np.round(mask_radii_res, 1)}
layer = dr+'-'+survey
#'dr8-model'
#'dr8-resid'
#layer2=layer+'-model'

masks = {'rmag<19':(BS_diff) & (bgsbutBS) & (cat['RMAG'] < 19), 'rmag>19':(BS_diff) & (bgsbutBS) & (cat['RMAG'] > 19)}

for key, val in zip(masks.keys(), masks.values()):
    
    idx = list(np.where(val))[0]
    print('sample size: %i' %(len(idx)))

    #if we want to compare the same objects with another catalogue (i.e., DR7, DR8) select only
    #the right number of indexes as your grid to avoid the random selection in postages_circle.
    postages_circle(coord=[cat['RA'], cat['DEC']], centeridx=idx, veto=veto, info=info, scale=0.262, 
                scale_unit='pixscale', layer=layer, radius=4/3600., m=2, grid=[12,6], 
                    savefile='%s/%s_BSold_not_BSnew_bgs_%s' %(pathdir, layer, key), layer2=layer+'-resid', layer2Mode='merge', 
                        isLG=False, title=None, markers=False, colorkey=True)
    
    plt.show()

Progress...N/A%|                                                    |
sample size: 240

Progress...100%|####################################################|

Colour key:
	 BGS --> lime
	 Not BGS --> royalblue
	 other --> red

Progress...N/A%|                                                    |
sample size: 937

Progress...100%|####################################################|

Colour key:
	 BGS --> lime
	 Not BGS --> royalblue
	 other --> red
In [44]:
from io_ import get_msmask

#load mask sources objects from SWEEPS DR8
masksources = np.load('/global/cscratch1/sd/qmxp55/bgstargets_output/dr8/dr8_sweep_whole_maskbitsources.npy')
In [45]:
#get the medium bright stars
starsMS = get_msmask(masksources)

6584 nearby objects
12378075 Bright Stars
In [ ]:
hppix_starsMS = hp.ang2pix(nside,(90.-starsMS['DEC'])*np.pi/180.,starsMS['RA']*np.pi/180.,nest=True) # catalogue hp pixels array
if reg[:8] == 'svfields': starsMSinreg = get_svfields_fg(starsMS['RA'], starsMS['DEC'])
else: starsMSinreg = get_reg(reg=reg, hppix=hppix_starsMS)

#cstars = SkyCoord(stars['RA']*units.degree,stars['DEC']*units.degree, frame='icrs')
#galb_stars = cstars.galactic.b.value # galb coordinate
starsMS = starsMS[(starsMSinreg) & (starsMS['DEC'] > -25)]
In [ ]:
#
log = True
nbins = 150

bgs_any = ((cat['BGSBITS'] & 2**(20)) != 0)
        
plt.figure(figsize=(20, 20))
plt.title(r'%s %s' %(dr, survey), size=20)
plt.axis('off')
_ = overdensity(cat[bgs_any], starsMS, BS_radii, 'MAG', 35, density=False, 
                        magbins=[13,16,3], radii_2=new_BS_radii, grid=[2,2], SR=[0., 110.], scaling=False, nbins=nbins, 
                               SR_scaling=3,logDenRat=[-5, 5], radii_bestfit=False, annulus=None, bintype='2', 
                                   filename='%s/2d_stack_MS_%s_%s' %(pathdir, dr, survey), log=log)

mag_bins_len: 3
13.00 < MAG < 14.00
mag_mean 13.554746475769027
mask_radius2 9.016305834215077
mag_radii MAX: 21.356534329212277 mag_radii MIN: 15.745354637543699
mag MAX: 13.999975204467773 mag MIN: 13.000082015991211
density cumu (min, max): (nan, nan)
density non-cumu (min, max): (-inf, 3.37)
2374 of inf in density ratio out of a total of 16977
Minimum density ratio = -5, Maximum density ratio = 2.42612
----------------
14.00 < MAG < 15.00
mag_mean 14.552994305402924
mask_radius2 6.650310656744731
mag_radii MAX: 15.744916668273756 mag_radii MIN: 11.600061258980851
mag MAX: 14.999937057495117 mag MIN: 14.000068664550781
density cumu (min, max): (nan, nan)
density non-cumu (min, max): (-inf, 4.07)
397 of inf in density ratio out of a total of 16977
Minimum density ratio = -5, Maximum density ratio = 1.84924
----------------
15.00 < MAG < 16.00
mag_mean 15.54923895213025
mask_radius2 4.900082750800707
mag_radii MAX: 11.599695559901352 mag_radii MIN: 8.544621042989482
mag MAX: 15.999944686889648 mag MIN: 15.00003719329834
density cumu (min, max): (nan, nan)
density non-cumu (min, max): (-inf, 2.69)
29 of inf in density ratio out of a total of 16977
Minimum density ratio = -5, Maximum density ratio = 1.75134
----------------

Grr vs g-z: bgsbutSG

In [ ]:
#bgslist = ['BS', 'LG', 'GC', 'nobs', 'FMC2', 'CC', 'QC_FM', 'QC_FI', 'QC_FF', 'QC_IVAR']
    
## get BGS objects without BS bit set on
#bgsbutSG = np.ones_like(cat['RA'], dtype=bool)

#for key in bgslist:
#    keep = ((cat['BGSBITS'] & 2**(bgsmask()[key])) != 0)
#    bgsbutSG &= keep
#bgsbutSG &= cat['RMAG'] < 20

bgsbutSG = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['SG'], bgsmask=bgsmask(), rlimit=20)

finite = np.ones_like(cat['RMAG'], dtype='?')
for i in ['RMAG', 'GMAG', 'ZMAG','FLUX_R', 'RFIBERMAG']:
    finite &= np.isfinite(cat[i])
In [ ]:
Ared = hpdict['bgsarea_'+reg]
In [ ]:
#
Grr = cat['G'] - flux_to_mag(cat['FLUX_R'])
coord = {'g-z':cat['GMAG'] - cat['ZMAG'], 'G-rr':Grr}
PSF = (cat['TYPE'] == 'PSF ') | (cat['TYPE'] == 'PSF') 
#if dr[:3] == 'dr9': morphos = ['REX', 'EXP', 'DEV', 'SER']
#else: morphos = ['REX ', 'EXP ', 'DEV ', 'COMP']
rows, cols = 2, 2
#debug = cat['RMAG'] < 16
mask = (catinreg) & (bgsbutSG) & (finite) #& (debug)
hline, vline = 0.6, None
vmin, vmax = 1, None
    
fig    = plt.figure(figsize=(8*cols,6*rows))
gs     = gridspec.GridSpec(rows, cols, hspace=0.2, wspace=0.10)
cmap = plt.get_cmap('coolwarm', 8)
    
for i, morpho in enumerate(morphos):
        
    if (i%cols==0): ylab=True
    else: ylab = False
    if (morpho == PSF): 
        vmax, cbar = None, 'vertical'
    else:
        vmax, cbar = 100, 'panel'
        
    if i < 2: xlab = False
    else: xlab = True
        
    morphomask = cat['TYPE'] == morpho
    hexbin(coord=coord, catmask=((mask) & (morphomask)), n=i, bins='log', title=morpho, cmap=cmap, 
               ylab=ylab, xlab=xlab, vline=vline, hline=hline, fig=fig, gs=gs, vmin=vmin, vmax=vmax, mincnt=1, 
                   file='%s/gz_Grr_bgsbutSG_hexbin_extended_%s' %(pathdir, survey), fracs=True, area=Ared, cbar=cbar)

fig    = plt.figure(figsize=(8*1,8*1))
gs     = gridspec.GridSpec(1, 1, hspace=0.15, wspace=0.10)

hexbin(coord=coord, catmask=((mask) & (PSF)), n=0, bins='log', title='PSF', cmap=cmap, 
               ylab=True, vline=vline, hline=hline, fig=fig, gs=gs, vmin=None, vmax=1000, mincnt=1, 
                   file='%s/gz_Grr_bgsbutSG_hexbin_psf_%s' %(pathdir, survey), fracs=True, area=Ared, cbar='horizontal')
Out[ ]:
<matplotlib.axes._subplots.AxesSubplot at 0x2aaae6b158d0>

rmag vs rfibmag: bgsbutFMC

In [ ]:
#bgslist = ['BS', 'LG', 'GC', 'nobs', 'SG', 'CC', 'QC_FM', 'QC_FI', 'QC_FF', 'QC_IVAR']
    
## get BGS objects without BS bit set on
#bgsbutFMC2 = np.ones_like(cat['RA'], dtype=bool)

#for key in bgslist:
#    keep = ((cat['BGSBITS'] & 2**(bgsmask()[key])) != 0)
#    bgsbutFMC2 &= keep
#bgsbutFMC2 &= cat['RMAG'] < 20

bgsbutFMC2 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['FMC2'], bgsmask=bgsmask(), rlimit=20)
In [ ]:
#
rows, cols = 2, 2
coord = {'r':cat['RMAG'], 'rfibmag':cat['RFIBERMAG']}

#debug = cat['RMAG'] < 16
mask = (catinreg) & (bgsbutFMC2) & (finite) #& (debug)
hline, vline = None, None
vmin, vmax = 1, 200
    
fig    = plt.figure(figsize=(8*cols,6*rows))
gs     = gridspec.GridSpec(rows, cols, hspace=0.25, wspace=0.10)
cmap = plt.get_cmap('coolwarm', 8)
    
for i, morpho in enumerate(morphos):
        
    if (i%cols==0): ylab=True
    else: ylab = False
    if morpho == PSF: 
        vmax, cbar = None, 'vertical'
    else:
        vmax, cbar = 2000, 'panel'
        
    if i < 2: xlab = False
    else: xlab = True
        
    morphomask = cat['TYPE'] == morpho
    hexbin(coord=coord, catmask=((mask) & (morphomask)), n=i, bins='log', title=morpho, cmap=cmap, 
               ylab=ylab, xlab=xlab, vline=vline, hline=hline, fig=fig, gs=gs, vmin=vmin, vmax=vmax, mincnt=1, fmcline=True, 
                   file='%s/r_rfibmag_bgsbutFMC2_hexbin_extended_%s' %(pathdir, survey), fracs=False, area=Ared, cbar=cbar)
    
    
fig    = plt.figure(figsize=(8*1,8*1))
gs     = gridspec.GridSpec(1, 1, hspace=0.15, wspace=0.10)

hexbin(coord=coord, catmask=((mask) & (PSF)), n=0, bins='log', title='PSF', cmap=cmap, 
               ylab=True, vline=vline, hline=hline, fig=fig, gs=gs, vmin=None, vmax=1000, mincnt=1, fmcline=True,
                   file='%s/r_rfibmag_bgsbutFMC2_hexbin_psf_%s' %(pathdir, survey), fracs=False, area=Ared, cbar='horizontal')
Out[ ]:
<matplotlib.axes._subplots.AxesSubplot at 0x2aaae577ac50>

g-r vs r-z: bgsbutCC

In [ ]:
bgsbutCC = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['CC'], bgsmask=bgsmask(), rlimit=20)
In [ ]:
#
import matplotlib.patches as patches

coord = {'g-r':cat['GMAG'] - cat['RMAG'], 'r-z':cat['RMAG'] - cat['ZMAG']}
mask = (catinreg) & (bgsbutCC) & (finite) # & (cat['RMAG'] < 14)

fig   = plt.figure(figsize=(10,10))
gs     = gridspec.GridSpec(1, 1, hspace=0.15, wspace=0.10)


hexbin(coord=coord, catmask=mask, n=0, bins='log', title=None, cmap=cmap, 
               ylab=True, vline=None, hline=None, fig=fig, gs=gs, xlim=(-2, 7), ylim=(-4, 5), 
                   vmin=None, vmax=None, mincnt=1, fmcline=False,
                       file=None, fracs=False, area=Ared, cbar='horizontal')

plt.plot(np.linspace(-1, 4, 3), np.full(3, -1), lw=2, c='r')
plt.plot(np.linspace(-1, 4, 3), np.full(3, 4), lw=2, c='r')

plt.plot(np.full(3, -1), np.linspace(-1, 4, 3), lw=2, c='r')
plt.plot(np.full(3, 4), np.linspace(-1, 4, 3), lw=2, c='r')

file='%s/gr_rz_bgsbutCC_hexbin_%s' %(pathdir, survey)
fig.savefig(file, bbox_inches = 'tight', pad_inches = 0)

Venn diagrams: QCs

In [ ]:
from QA import plot_venn3
In [ ]:
# get BGS but QCs
bgsbutQCs = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['QC_FM', 'QC_FI', 'QC_FF'], bgsmask=bgsmask(), rlimit=20)
FM = ((cat['BGSBITS'] & 2**(11)) == 0) #rejects by FM
FI = ((cat['BGSBITS'] & 2**(12)) == 0) #rejects by FI
FF = ((cat['BGSBITS'] & 2**(13)) == 0) #rejects by FF
area = hpdict0['bgsarea_'+survey]
In [ ]:
#
keepB = (bgsbutQCs) & (cat['RMAG'] < 19.5)
keepF = (bgsbutQCs) & (cat['RMAG'] > 19.5) & (cat['RMAG'] < 20)
plot_venn3(A=(keepB & FM), B=(keepB & FI), C=(keepB & FF), norm=area, 
           labels=['FM', 'FI', 'FF'], file=pathdir+'/venn_QCs_bright_%s' %(survey), title='BRIGHT', colors = ['deepskyblue', 'springgreen', 'salmon'])
plot_venn3(A=(keepF & FF), B=(keepF & FI), C=(keepF & FM), norm=area, 
           labels=['FF', 'FI', 'FM'], file=pathdir+'/venn_QCs_faint_%s' %(survey), title='FAINT', colors = ['salmon' ,'springgreen', 'deepskyblue'])

AEN vs G-rr star-galaxy classification in GAIA

Tractor force to fit as PSF the GAIA objects that satisfy the following conditions:

  • For $G<19$: astrometric_excess_noise $< 10^{0.5}$,

  • For $G≥19$: astrometric_excess_noise $ < 10^{0.5+0.2(G−19)}$.

G is the GAIA photometric $G$-band and the astrometric excess noise ($AEN$) is the... Those PSF that are in GAIA are treated as Stars according to TRACTOR and the GAIA objects that do not fulfill that condition are treated as GALAXIES. This is different to the way we do star-galaxy separation for BGS, we use GAIA but rather than look at the $AEN$ we look at the colour distribution $G-rr$.

The PSF GAIA objects (from above definition) are treated differently than other PSF objects and they have a local fitting, the problem is that maybe those are not stars and might be galaxies instead that photometry is being altered by this local fitting.

Here we are going to find out what are the differences between our star-galaxy classification and the one above for GAIA objects. We're going to:

  • $Grr$ distributions for $Grr$ class and $AEN$ class:
    • All
    • TRACTOR PSF
    • TRACTOR No PSF
  • Heat maps
  • Venn diagram's
  • Galleries of mismatches divided in TRACTOR PSF and TRACTOR non-PSF:
    • Stars with $Grr$ class and Galaxies with $AEN$ class
    • Galaxies with $Grr$ class and Stars with $AEN$ class
In [ ]:
inGAIA = (cat['REF_CAT'] == 'G2') & (cat['RMAG'] < 20)
grr_gal = ((cat['BGSBITS'] & 2**(6)) != 0) & (cat['RMAG'] < 20)
grr_stars = ((cat['BGSBITS'] & 2**(6)) == 0) & (cat['RMAG'] < 20)
AEN_star, AEN_gal = gaiaAEN(inGAIA=inGAIA, size=len(cat['RA']), G=cat['G'], AEN=cat['AEN'], dr=dr)
PSF = (cat['TYPE'] == 'PSF ') | (cat['TYPE'] == 'PSF')

#inGAIA2 = (cat['G'] != 0)
In [ ]:
#
coord = {'G':cat['G'], r'$\log(AEN)$':np.log10(cat['AEN'])}
#masks = [(grr_stars) & (inGAIA) & (cat['RMAG'] < 14), (grr_gal) & (inGAIA) & (cat['RMAG'] < 14)]
masks = [(grr_stars) & (inGAIA), (grr_gal) & (inGAIA)]
titles = ['stars with G-rr (r < 20)', 'galaxies with G-rr (r < 20)']
vmin, vmax = 1, None
    
fig    = plt.figure(figsize=(9*len(masks), 7))
gs     = gridspec.GridSpec(1, len(masks),hspace=0.1,wspace=0.10)
cmap = plt.get_cmap('coolwarm', 8)
xlim=[10, 22]
ylim=[-2, 1.5]

for i, mask in enumerate(masks):
        
    #if (i%len(masks)==0): xlab=True
    #else: xlab = False
    #if morpho == 'PSF ': vmax = None
        
    #morphomask = cat['TYPE'] == morpho
    hexbin(coord=coord, catmask=mask, n=i, bins='log', title=titles[i], cmap=cmap, xlab=True,
               ylab=True, vline=None, hline=None, fig=fig, gs=gs, vmin=vmin, vmax=vmax, mincnt=1, fmcline=False, 
                   file=None, area=None, xlim=xlim, ylim=ylim, cbar='horizontal')
    
    #if i == 0:
    #    x, y = coord.keys()
    #    mask1 = (gaia_star) & (dr8_gama) & (dr8_z > 0.002) & (bgs)
    #    mask2 = (gaia_star) & (dr8_gama) & (dr8_z > 0.002) & (~bgs)
    #    plt.scatter(coord[x][mask2], coord[y][mask2], s=2, c='royalblue', label=r'Not in BGS')
    #    plt.scatter(coord[x][mask1], coord[y][mask1], s=2, c='k', label=r'in BGS')
    #    plt.legend()
    
    x_ = np.linspace(19, xlim[1], 4)
    plt.plot(np.linspace(xlim[0], 18, 3), np.full(3, 0.5), color='r', ls='-', lw=2)
    plt.plot(np.full(3, 18), np.linspace(ylim[0], 0.5, 3), color='r', ls='-', lw=2)
    plt.plot(x_, 0.5 + 0.2*(x_ - 19.), color='r', ls='--', lw=2, alpha=0.7)
    plt.plot(np.linspace(18,19,3), np.full(3, 0.5), color='r', ls='--', lw=2, alpha=0.7)

file = pathdir+'/G_logAEN'
fig.savefig(file+'.png', bbox_inches = 'tight', pad_inches = 0)
  • How many galaxies (i.e. G-rr > 0.6) that are AEN stars? Are these real galaxies?
  • How many SDSS & GAMA "truth" galaxies are AEN star? plot SDSS & GAMA galaxies on top of above plot.
    • how many of above in G-rr > 0.6?
    • how many of above in G-rr < 0.6?
In [ ]:
#
area = hpdict0['bgsarea_'+survey]
BS = ((cat['BGSBITS'] & 2**(0)) != 0)
gal_in_aenstar = (inGAIA) & (grr_gal) & (cat['G'] < 18)*(np.log10(cat['AEN']) < 0.5)
print('G-rr > 0.6 & in GAIA & G < 18 & AEN < 10**(0.5): \t %i, \t %.3f (1/deg^2)' 
      %(np.sum(gal_in_aenstar), np.sum(gal_in_aenstar)/area))
print('G-rr > 0.6 & in GAIA & G < 18 & AEN < 10**(0.5) & BS: \t %i, \t %.3f (1/deg^2)' 
      %(np.sum(gal_in_aenstar & (BS)), np.sum(gal_in_aenstar & (BS))/area))
print('G-rr > 0.6 & in GAIA & G < 18 & AEN < 10**(0.5) & BGS: \t %i, \t %.3f (1/deg^2)' 
      %(np.sum(gal_in_aenstar & (bgs_any)), np.sum(gal_in_aenstar & (bgs_any))/area))

G-rr > 0.6 & in GAIA & G < 18 & AEN < 10**(0.5): 	 440, 	 3.962 (1/deg^2)
G-rr > 0.6 & in GAIA & G < 18 & AEN < 10**(0.5) & BS: 	 158, 	 1.423 (1/deg^2)
G-rr > 0.6 & in GAIA & G < 18 & AEN < 10**(0.5) & BGS: 	 12, 	 0.108 (1/deg^2)
In [61]:
#
veto = {'BGS':(bgs_any),'Not BGS':(~bgs_any)}

info = {'r':cat['RMAG'], 'T':cat['TYPE'].astype(str), 'ref':cat['REF_CAT'].astype(str)}
layer = dr+'-'+survey
#'dr8-model'
#'dr8-resid'
#layer2=layer+'-model'

idx = list(np.where((gal_in_aenstar) & (bgs_any)))[0]
print('sample size: %i' %(len(idx)))

#if we want to compare the same objects with another catalogue (i.e., DR7, DR8) select only
#the right number of indexes as your grid to avoid the random selection in postages_circle.
postages_circle(coord=[cat['RA'], cat['DEC']], centeridx=idx, veto=veto, info=info, scale=0.262, 
                scale_unit='pixscale', layer=layer, radius=4/3600., m=2, grid=[2,5], 
                    savefile='%s/%s_Grr_gal_AEN_stars' %(pathdir, layer), layer2=layer+'-resid', layer2Mode='merge', 
                        isLG=False, title=None, markers=False, colorkey=True)

Progress...N/A%|                                                    |
sample size: 12

Progress...100%|####################################################|

Colour key:
	 BGS --> lime
	 Not BGS --> royalblue
	 other --> red
In [62]:
#
fig = plt.figure(figsize=(18,12))

log = True
Grr = cat['G'] - flux_to_mag(cat['FLUX_R'])

classes = {'stars_Grr':(grr_stars) & (inGAIA), 'gal_Grr':(grr_gal) & (inGAIA), 'stars_aen': (AEN_star) & (cat['RMAG'] < 20), 'gal_aen':(AEN_gal) & (cat['RMAG'] < 20), 
              'stars_Grr_PSF':(grr_stars) & (inGAIA) & (PSF), 'gal_Grr_PSF':(grr_gal) & (inGAIA) & (PSF),
                  'stars_aen_PSF':(AEN_star) & (PSF) & (cat['RMAG'] < 20), 'gal_aen_PSF':(AEN_gal) & (PSF) & (cat['RMAG'] < 20),
                      'stars_Grr_noPSF':(grr_stars) & (inGAIA) & (~PSF), 'gal_Grr_noPSF':(grr_gal) & (inGAIA) & (~PSF),
                          'stars_aen_noPSF':(AEN_star) & (~PSF) & (cat['RMAG'] < 20), 'gal_aen_noPSF':(AEN_gal) & (~PSF) & (cat['RMAG'] < 20)}

for i, j in enumerate(['Grr', 'aen', 'Grr_PSF', 'aen_PSF', 'Grr_noPSF', 'aen_noPSF']):

    XS = (classes['stars_' + j])
    XG = (classes['gal_' + j])
    tot_gaia = np.sum((inGAIA))

    plt.subplot(3, 2, i+1)
    
    if i == 0:
        plt.title(r'$G-rr$ classification (r < 20)', size=22)
    if i == 1:
        plt.title(r'$AEN$ classification (r < 20)', size=22)
    
    bins = np.linspace(-2.1, 5, 60)
        
    plt.hist(Grr[XS], bins=bins, log=log, alpha=0.7, label='%s stars, $f$=%2.3g' %(j, np.sum(XS)/tot_gaia))
    plt.hist(Grr[XG], bins=bins, log=log, alpha=0.4, label='%s galaxies, $f$=%2.3g' %(j, np.sum(XG)/tot_gaia))
    if i == 0:
        plt.axvline(0.6, ls='--', lw=2, c='k', label='SG threshold, Grr = 0.6')
        plt.text(-2.1, 2*10**7, '$N_{tot}$=%s' %(tot_gaia))
    plt.axvline(0.6, ls='--', lw=2, c='k')
    if (i == 4) or (i == 5):
        plt.xlabel(r'$G-rr$', size=20)
    if i == 2:
        plt.ylabel(r'$N$', size=20)    
    plt.legend()

plt.show()
filename = '%s/gaia_Grr_comparison_with_aen' %(pathdir)
fig.savefig(filename+'.png') #bbox_inches = 'tight', pad_inches = 0
In [63]:
from matplotlib_venn import venn2, venn2_circles

fig = plt.figure(figsize=(5,5))

sf = 2
area = hpdict0['bgsarea_'+survey]

a = classes['gal_Grr']
b = classes['gal_aen']
c = (a) & (b)

a1 = round(((np.sum(a) - np.sum(c))/area), sf)
b1 = round(((np.sum(b) - np.sum(c))/area), sf)
c1 = round(np.sum(c)/area, sf)

plt.title(r'GALAXIES', size=18)
labels = ('Grr', 'AEN')
   
venn2([a1, b1, c1], set_labels = labels)
c=venn2_circles([a1, b1, c1], linestyle='solid', linewidth=1, color="k")

filename = '%s/venn_gaia_Grr_comparison_with_aen' %(pathdir)
fig.savefig(filename+'.png', bbox_inches = 'tight', pad_inches = 0)

Systematics

In [59]:
#
if reg[:8] == 'svfields': reg_ = reg+'_'+survey[0]
else: reg_ = reg

fig    = plt.figure(figsize=(18,5))
gs     = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)
isdesi     = (hpdict['isdesi']) & (hpdict['bgsfracarea']>0) #this should be rename as "mainreg"
nx         = 20

b0, m0 = plot_sysdens(hpdicttmp=hpdict, namesels=['bgs_any'], regs=[reg_], syst='stardens', mainreg=isdesi, 
                      xlim=None, n=0, nx=nx, clip=True, denslims=False, ylab=True, weights=True, 
                          fig=fig, gs=gs, title='test', label=True)
In [49]:
#dic WEIGHTED

ws = 1./((m0)*hpdict['stardens'] + b0)

hpdict_ws = get_dict(cat=cat, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=hppix_cat, 
                      maskrand=((rancuts['BS']) & (rancuts['GC']) & (rancuts['LG']) & (rancuts['nobs'])),
                          maskcat=None, 
                             Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None, 
                                desifootprint=True, namesels={'bgs_any':20, 'bright':21, 'faint':22}, 
                                  galb=cat_ex['b'], log=True, survey='bgs', ws=ws)

positions and desifotprint DONE...
systematics DONE...
randdens =  15000  ; len randoms =  68230483
bgsfracarea DONE...
regions DONE...
areas DONE...
computing for  north / bgs_any
computing for  north / bright
computing for  north / faint
target densities in north DONE...
computing for  south / bgs_any
computing for  south / bright
computing for  south / faint
target densities in south DONE...
meandens_bgs_any_all = 1398 /deg2
meandens_bgs_any_des = nan /deg2
meandens_bgs_any_decals = nan /deg2
meandens_bgs_any_north = 1398 /deg2
meandens_bgs_any_south = nan /deg2
meandens_bgs_any_svfields = 1370 /deg2
meandens_bgs_any_svfields_n = 1370 /deg2
meandens_bgs_any_svfields_s = nan /deg2
meandens_bright_all = 824 /deg2
meandens_bright_des = nan /deg2
meandens_bright_decals = nan /deg2
meandens_bright_north = 824 /deg2
meandens_bright_south = nan /deg2
meandens_bright_svfields = 798 /deg2
meandens_bright_svfields_n = 798 /deg2
meandens_bright_svfields_s = nan /deg2
meandens_faint_all = 574 /deg2
meandens_faint_des = nan /deg2
meandens_faint_decals = nan /deg2
meandens_faint_north = 574 /deg2
meandens_faint_south = nan /deg2
meandens_faint_svfields = 573 /deg2
meandens_faint_svfields_n = 573 /deg2
meandens_faint_svfields_s = nan /deg2
dens_bgs_any_all = 1397 /deg2
dens_bgs_any_des = nan /deg2
dens_bgs_any_decals = nan /deg2
dens_bgs_any_north = 1397 /deg2
dens_bgs_any_south = nan /deg2
dens_bgs_any_svfields = 1369 /deg2
dens_bgs_any_svfields_n = 1369 /deg2
dens_bgs_any_svfields_s = nan /deg2
dens_bright_all = 823 /deg2
dens_bright_des = nan /deg2
dens_bright_decals = nan /deg2
dens_bright_north = 823 /deg2
dens_bright_south = nan /deg2
dens_bright_svfields = 796 /deg2
dens_bright_svfields_n = 796 /deg2
dens_bright_svfields_s = nan /deg2
dens_faint_all = 574 /deg2
dens_faint_des = nan /deg2
dens_faint_decals = nan /deg2
dens_faint_north = 574 /deg2
dens_faint_south = nan /deg2
dens_faint_svfields = 572 /deg2
dens_faint_svfields_n = 572 /deg2
dens_faint_svfields_s = nan /deg2
In [ ]:
In [50]:
from QA import pixhistregs
In [51]:
# dr8_south+north : density distributions + systematics

# settings
isdesi     = (hpdict['isdesi']) & (hpdict['bgsfracarea']>0) #this should be rename as "mainreg"
systs      = ['log10_stardens','ebv','psfsize_g', 'psfsize_r', 'psfsize_z','galdepth_g', 'galdepth_r', 'galdepth_z']
nx         = 20
cols       = ['0.5','b','g','r']    
#regs = [reg]
namesels=['bgs_any', 'bright', 'faint']
#

# looping on subselections
#for key, title in zip([None, ws], ['UNWEIGHTED', 'WEIGHTED']):
for key, title in zip([hpdict, hpdict_ws], ['UNWEIGHTED', 'WEIGHTED']):
    
    # systematics
    fig    = plt.figure(figsize=(20,15))
    gs     = gridspec.GridSpec(3,3,hspace=0.30,wspace=0.25)
        
    for i in range(len(systs)+1):
        
        if i == 0:
            axinfo = fig.add_subplot(gs[i])
            # infos
            axinfo.axis('off')
            #axinfo.text(0.5,0.9,namesel,fontsize=20,fontweight='bold',ha='center',transform=axinfo.transAxes)
            tmpy = 0.7
            #for regi,col in zip(regs,cols):
            #if reg == 'south':reg_ = 'DECaLS+DES'
            #else: reg_ = reg
            for namesel,col in zip(namesels,cols):
                tmpstr = surveylab+' & '+namesel+' : '+r'$\overline{n}=$'+'%.0f'%key['meandens_'+namesel+'_'+reg]+r' deg$^{-2}$'
                axinfo.text(0.05,tmpy,tmpstr,color=col,fontsize=15,transform=axinfo.transAxes)
                tmpy  -= 0.15
        else:
            syst = systs[i-1]
            if (i%3==0) or (i==1): ylab=True
            else: ylab = False
            if i == 1: label = True
            else: label = False
            #if title == 'UNWEIGHTED': 
            #    weights = False
            #    onlyweights=False
            #else: 
            #    weights = True
            #    onlyweights=True
                
            plot_sysdens(hpdicttmp=key, namesels=namesels, regs=[reg], syst=syst, mainreg=isdesi, xlim=None, n=i, nx=nx, clip=True, 
                         denslims=False, ylab=ylab, weights=False, fig=fig, gs=gs, label=label, title=title, 
                             ws=None)
    
    # save fig
    fig.savefig('%s/systematics_main_bgs_%s_%s.png' %(pathdir, survey, title), bbox_inches = 'tight', pad_inches = 0)
    print('')
    print('')

Systematics separated by best model fitted

In [19]:
namesels = {}
for i in ['COMP', 'DEV ', 'EXP ', 'PSF ', 'REX ']:
    namesels[i] = (cat['TYPE'] == i) & (bgs_any)
In [20]:
hpdict_tmp = get_dict(cat=cat, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=hppix_cat, 
                      maskrand=((rancuts['BS']) & (rancuts['GC']) & (rancuts['LG']) & (rancuts['nobs'])),
                          maskcat=None, 
                             Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None, 
                                desifootprint=True, namesels=namesels, galb=cat_ex['b'], log=True, survey='custom')

positions and desifotprint DONE...
systematics DONE...
randdens =  15000  ; len randoms =  145754705
bgsfracarea DONE...
regions DONE...
areas DONE...
computing for  north / COMP
computing for  north / DEV 
computing for  north / EXP 
computing for  north / PSF 
computing for  north / REX 
target densities in north DONE...
computing for  south / COMP
computing for  south / DEV 
computing for  south / EXP 
computing for  south / PSF 
computing for  south / REX 
target densities in south DONE...
meandens_COMP_all = 31 /deg2
meandens_COMP_des = 50 /deg2
meandens_COMP_decals = 28 /deg2
meandens_COMP_north = nan /deg2
meandens_COMP_south = 31 /deg2
meandens_COMP_svfields = 34 /deg2
meandens_COMP_svfields_n = nan /deg2
meandens_COMP_svfields_s = 34 /deg2
meandens_DEV _all = 630 /deg2
meandens_DEV _des = 650 /deg2
meandens_DEV _decals = 627 /deg2
meandens_DEV _north = nan /deg2
meandens_DEV _south = 630 /deg2
meandens_DEV _svfields = 655 /deg2
meandens_DEV _svfields_n = nan /deg2
meandens_DEV _svfields_s = 655 /deg2
meandens_EXP _all = 514 /deg2
meandens_EXP _des = 526 /deg2
meandens_EXP _decals = 513 /deg2
meandens_EXP _north = nan /deg2
meandens_EXP _south = 514 /deg2
meandens_EXP _svfields = 520 /deg2
meandens_EXP _svfields_n = nan /deg2
meandens_EXP _svfields_s = 520 /deg2
meandens_PSF _all = 5 /deg2
meandens_PSF _des = 2 /deg2
meandens_PSF _decals = 5 /deg2
meandens_PSF _north = nan /deg2
meandens_PSF _south = 5 /deg2
meandens_PSF _svfields = 4 /deg2
meandens_PSF _svfields_n = nan /deg2
meandens_PSF _svfields_s = 4 /deg2
meandens_REX _all = 246 /deg2
meandens_REX _des = 201 /deg2
meandens_REX _decals = 252 /deg2
meandens_REX _north = nan /deg2
meandens_REX _south = 246 /deg2
meandens_REX _svfields = 241 /deg2
meandens_REX _svfields_n = nan /deg2
meandens_REX _svfields_s = 241 /deg2
dens_COMP_all = 31 /deg2
dens_COMP_des = 50 /deg2
dens_COMP_decals = 28 /deg2
dens_COMP_north = nan /deg2
dens_COMP_south = 31 /deg2
dens_COMP_svfields = 34 /deg2
dens_COMP_svfields_n = nan /deg2
dens_COMP_svfields_s = 34 /deg2
dens_DEV _all = 629 /deg2
dens_DEV _des = 649 /deg2
dens_DEV _decals = 626 /deg2
dens_DEV _north = nan /deg2
dens_DEV _south = 629 /deg2
dens_DEV _svfields = 654 /deg2
dens_DEV _svfields_n = nan /deg2
dens_DEV _svfields_s = 654 /deg2
dens_EXP _all = 514 /deg2
dens_EXP _des = 525 /deg2
dens_EXP _decals = 512 /deg2
dens_EXP _north = nan /deg2
dens_EXP _south = 514 /deg2
dens_EXP _svfields = 519 /deg2
dens_EXP _svfields_n = nan /deg2
dens_EXP _svfields_s = 519 /deg2
dens_PSF _all = 5 /deg2
dens_PSF _des = 2 /deg2
dens_PSF _decals = 5 /deg2
dens_PSF _north = nan /deg2
dens_PSF _south = 5 /deg2
dens_PSF _svfields = 4 /deg2
dens_PSF _svfields_n = nan /deg2
dens_PSF _svfields_s = 4 /deg2
dens_REX _all = 246 /deg2
dens_REX _des = 201 /deg2
dens_REX _decals = 252 /deg2
dens_REX _north = nan /deg2
dens_REX _south = 246 /deg2
dens_REX _svfields = 240 /deg2
dens_REX _svfields_n = nan /deg2
dens_REX _svfields_s = 240 /deg2
In [27]:
fig    = plt.figure(figsize=(18,5))
gs     = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)
isdesi     = (hpdict_tmp['isdesi']) & (hpdict_tmp['bgsfracarea']>0)

plot_sysdens(hpdicttmp=hpdict_tmp, namesels=['COMP', 'DEV ', 'EXP ', 'REX '], regs=[reg], 
                 syst='log10_stardens', mainreg=isdesi, xlim=None, n=0, nx=20, clip=True, 
                         denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test2', label=True)

Systematics separated by rmag bins

In [30]:
namesels = {}
mags = np.linspace(16, 20, 9)
for i in range(len(mags[:-1])):
    namesels['%s_%s' %(str(mags[i]), str(mags[i+1]))] = (cat['RMAG'] > mags[i]) & (cat['RMAG'] < mags[i+1]) & (bgs_any)
In [31]:
namesels.keys()
Out[31]:
dict_keys(['16.0_16.5', '16.5_17.0', '17.0_17.5', '17.5_18.0', '18.0_18.5', '18.5_19.0', '19.0_19.5', '19.5_20.0'])
In [32]:
hpdict_tmp = get_dict(cat=cat, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=hppix_cat, 
                      maskrand=((rancuts['BS']) & (rancuts['GC']) & (rancuts['LG']) & (rancuts['nobs'])),
                          maskcat=None, 
                             Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None, 
                                desifootprint=True, namesels=namesels, galb=cat_ex['b'], log=True, survey='custom')

positions and desifotprint DONE...
systematics DONE...
randdens =  15000  ; len randoms =  145754705
bgsfracarea DONE...
regions DONE...
areas DONE...
computing for  north / 16.0_16.5
computing for  north / 16.5_17.0
computing for  north / 17.0_17.5
computing for  north / 17.5_18.0
computing for  north / 18.0_18.5
computing for  north / 18.5_19.0
computing for  north / 19.0_19.5
computing for  north / 19.5_20.0
target densities in north DONE...
computing for  south / 16.0_16.5
computing for  south / 16.5_17.0
computing for  south / 17.0_17.5
computing for  south / 17.5_18.0
computing for  south / 18.0_18.5
computing for  south / 18.5_19.0
computing for  south / 19.0_19.5
computing for  south / 19.5_20.0
target densities in south DONE...
meandens_16.0_16.5_all = 12 /deg2
meandens_16.0_16.5_des = 12 /deg2
meandens_16.0_16.5_decals = 12 /deg2
meandens_16.0_16.5_north = nan /deg2
meandens_16.0_16.5_south = 12 /deg2
meandens_16.0_16.5_svfields = 11 /deg2
meandens_16.0_16.5_svfields_n = nan /deg2
meandens_16.0_16.5_svfields_s = 11 /deg2
meandens_16.5_17.0_all = 21 /deg2
meandens_16.5_17.0_des = 21 /deg2
meandens_16.5_17.0_decals = 21 /deg2
meandens_16.5_17.0_north = nan /deg2
meandens_16.5_17.0_south = 21 /deg2
meandens_16.5_17.0_svfields = 20 /deg2
meandens_16.5_17.0_svfields_n = nan /deg2
meandens_16.5_17.0_svfields_s = 20 /deg2
meandens_17.0_17.5_all = 39 /deg2
meandens_17.0_17.5_des = 39 /deg2
meandens_17.0_17.5_decals = 39 /deg2
meandens_17.0_17.5_north = nan /deg2
meandens_17.0_17.5_south = 39 /deg2
meandens_17.0_17.5_svfields = 39 /deg2
meandens_17.0_17.5_svfields_n = nan /deg2
meandens_17.0_17.5_svfields_s = 39 /deg2
meandens_17.5_18.0_all = 71 /deg2
meandens_17.5_18.0_des = 71 /deg2
meandens_17.5_18.0_decals = 71 /deg2
meandens_17.5_18.0_north = nan /deg2
meandens_17.5_18.0_south = 71 /deg2
meandens_17.5_18.0_svfields = 70 /deg2
meandens_17.5_18.0_svfields_n = nan /deg2
meandens_17.5_18.0_svfields_s = 70 /deg2
meandens_18.0_18.5_all = 124 /deg2
meandens_18.0_18.5_des = 124 /deg2
meandens_18.0_18.5_decals = 124 /deg2
meandens_18.0_18.5_north = nan /deg2
meandens_18.0_18.5_south = 124 /deg2
meandens_18.0_18.5_svfields = 124 /deg2
meandens_18.0_18.5_svfields_n = nan /deg2
meandens_18.0_18.5_svfields_s = 124 /deg2
meandens_18.5_19.0_all = 211 /deg2
meandens_18.5_19.0_des = 213 /deg2
meandens_18.5_19.0_decals = 211 /deg2
meandens_18.5_19.0_north = nan /deg2
meandens_18.5_19.0_south = 211 /deg2
meandens_18.5_19.0_svfields = 217 /deg2
meandens_18.5_19.0_svfields_n = nan /deg2
meandens_18.5_19.0_svfields_s = 217 /deg2
meandens_19.0_19.5_all = 354 /deg2
meandens_19.0_19.5_des = 356 /deg2
meandens_19.0_19.5_decals = 353 /deg2
meandens_19.0_19.5_north = nan /deg2
meandens_19.0_19.5_south = 354 /deg2
meandens_19.0_19.5_svfields = 364 /deg2
meandens_19.0_19.5_svfields_n = nan /deg2
meandens_19.0_19.5_svfields_s = 364 /deg2
meandens_19.5_20.0_all = 579 /deg2
meandens_19.5_20.0_des = 580 /deg2
meandens_19.5_20.0_decals = 578 /deg2
meandens_19.5_20.0_north = nan /deg2
meandens_19.5_20.0_south = 579 /deg2
meandens_19.5_20.0_svfields = 596 /deg2
meandens_19.5_20.0_svfields_n = nan /deg2
meandens_19.5_20.0_svfields_s = 596 /deg2
dens_16.0_16.5_all = 12 /deg2
dens_16.0_16.5_des = 12 /deg2
dens_16.0_16.5_decals = 12 /deg2
dens_16.0_16.5_north = nan /deg2
dens_16.0_16.5_south = 12 /deg2
dens_16.0_16.5_svfields = 11 /deg2
dens_16.0_16.5_svfields_n = nan /deg2
dens_16.0_16.5_svfields_s = 11 /deg2
dens_16.5_17.0_all = 21 /deg2
dens_16.5_17.0_des = 21 /deg2
dens_16.5_17.0_decals = 21 /deg2
dens_16.5_17.0_north = nan /deg2
dens_16.5_17.0_south = 21 /deg2
dens_16.5_17.0_svfields = 20 /deg2
dens_16.5_17.0_svfields_n = nan /deg2
dens_16.5_17.0_svfields_s = 20 /deg2
dens_17.0_17.5_all = 39 /deg2
dens_17.0_17.5_des = 39 /deg2
dens_17.0_17.5_decals = 39 /deg2
dens_17.0_17.5_north = nan /deg2
dens_17.0_17.5_south = 39 /deg2
dens_17.0_17.5_svfields = 39 /deg2
dens_17.0_17.5_svfields_n = nan /deg2
dens_17.0_17.5_svfields_s = 39 /deg2
dens_17.5_18.0_all = 71 /deg2
dens_17.5_18.0_des = 71 /deg2
dens_17.5_18.0_decals = 71 /deg2
dens_17.5_18.0_north = nan /deg2
dens_17.5_18.0_south = 71 /deg2
dens_17.5_18.0_svfields = 70 /deg2
dens_17.5_18.0_svfields_n = nan /deg2
dens_17.5_18.0_svfields_s = 70 /deg2
dens_18.0_18.5_all = 124 /deg2
dens_18.0_18.5_des = 124 /deg2
dens_18.0_18.5_decals = 124 /deg2
dens_18.0_18.5_north = nan /deg2
dens_18.0_18.5_south = 124 /deg2
dens_18.0_18.5_svfields = 124 /deg2
dens_18.0_18.5_svfields_n = nan /deg2
dens_18.0_18.5_svfields_s = 124 /deg2
dens_18.5_19.0_all = 211 /deg2
dens_18.5_19.0_des = 213 /deg2
dens_18.5_19.0_decals = 211 /deg2
dens_18.5_19.0_north = nan /deg2
dens_18.5_19.0_south = 211 /deg2
dens_18.5_19.0_svfields = 216 /deg2
dens_18.5_19.0_svfields_n = nan /deg2
dens_18.5_19.0_svfields_s = 216 /deg2
dens_19.0_19.5_all = 353 /deg2
dens_19.0_19.5_des = 356 /deg2
dens_19.0_19.5_decals = 353 /deg2
dens_19.0_19.5_north = nan /deg2
dens_19.0_19.5_south = 353 /deg2
dens_19.0_19.5_svfields = 363 /deg2
dens_19.0_19.5_svfields_n = nan /deg2
dens_19.0_19.5_svfields_s = 363 /deg2
dens_19.5_20.0_all = 578 /deg2
dens_19.5_20.0_des = 579 /deg2
dens_19.5_20.0_decals = 578 /deg2
dens_19.5_20.0_north = nan /deg2
dens_19.5_20.0_south = 578 /deg2
dens_19.5_20.0_svfields = 595 /deg2
dens_19.5_20.0_svfields_n = nan /deg2
dens_19.5_20.0_svfields_s = 595 /deg2
In [33]:
fig    = plt.figure(figsize=(18,5))
gs     = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)
isdesi     = (hpdict_tmp['isdesi']) & (hpdict_tmp['bgsfracarea']>0)

plot_sysdens(hpdicttmp=hpdict_tmp, namesels=['16.0_16.5', '16.5_17.0', '17.0_17.5', '17.5_18.0'], regs=[reg], 
                 syst='log10_stardens', mainreg=isdesi, xlim=None, n=0, nx=20, clip=True, 
                         denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test3', label=True)
In [34]:
fig    = plt.figure(figsize=(18,5))
gs     = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)
isdesi     = (hpdict_tmp['isdesi']) & (hpdict_tmp['bgsfracarea']>0)

plot_sysdens(hpdicttmp=hpdict_tmp, namesels=['18.0_18.5', '18.5_19.0', '19.0_19.5', '19.5_20.0'], regs=[reg], 
                 syst='log10_stardens', mainreg=isdesi, xlim=None, n=0, nx=20, clip=True, 
                         denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test4', label=True)

Skymap footprint

Below code only needs the randoms.

In [5]:
#
from QA import set_mwd, get_radec_mw
from io_ import get_isdes

# dict without geometrical mask and in WHOLE SWEEPS footprint
hpdict0 = get_dict(cat=None, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=None, 
                      maskrand=None, maskcat=None, Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None, 
                            desifootprint=False, namesels=None, target_outputs=False, log=False)
# dict without geometrical mask and in DESI footprint
hpdict1 = get_dict(cat=None, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=None, 
                      maskrand=None, maskcat=None, Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None, 
                            desifootprint=True, namesels=None, target_outputs=False, log=False)

#
for i in ['area_all', 'bgsarea_south', 'bgsarea_north', 'bgsarea_des']:
    print('%s: \t %.3f(whole) \t %.3f(desifootprint)' %(i, hpdict0[i], hpdict1[i]))
    
if survey == 'south':
    # DECaLS area in the NGC/SGC within DESI footprint
    area_decals_ngc = hpdict1['bgsfracarea'][(hpdict1['issouth']) & (hpdict1['galb'] > 0)].sum() * pixarea
    area_decals_sgc = hpdict1['bgsfracarea'][(hpdict1['issouth']) & (hpdict1['galb'] < 0)].sum() * pixarea
    print('NGC DECaLS in DESI footprint: \t %.3f' %(area_decals_ngc))
    print('SGC DECaLS in DESI footprint: \t %.3f' %(area_decals_sgc))
    
if survey == 'north':
    # DECaLS area in the NGC/SGC within DESI footprint
    area_bass_ngc = hpdict1['bgsfracarea'][(hpdict1['isnorth']) & (hpdict1['galb'] > 0)].sum() * pixarea
    print('NGC BASS/MzLS in DESI footprint: \t %.3f' %(area_bass_ngc))

desitiles = fitsio.read('/global/homes/q/qmxp55/desi-tiles-viewer.fits')
org          = 120  # centre ra for mollweide plots
projection   = 'mollweide'

# plotting bgsfracarea
fig = plt.figure(figsize=(14, 7))
gs = gridspec.GridSpec(1,1,wspace=0.15,hspace=0)

ax = plt.subplot(gs[0],projection=projection)
_ = set_mwd(ax,org=org)
ramw,decmw = get_radec_mw(hpdict0['ra'],hpdict0['dec'],org)
tmp        = (hpdict0['bgsfracarea']>0) & (hpdict0['is'+survey])
if survey == 'south':
    SC         = ax.scatter(ramw[tmp],decmw[tmp],s=0.4, color='gray', label=r'DECam imaging')
if survey == 'north':
    SC         = ax.scatter(ramw[tmp],decmw[tmp],s=0.4, color='gray', label=r'BASS/MzLS imaging')

ramw,decmw = get_radec_mw(desitiles['ra'],desitiles['dec'],org)
if survey == 'south':
    tmp = (desitiles['in_desi'] ==  1) & (desitiles['dec'] < dec_resol_ns)
    SC         = ax.scatter(ramw[tmp],decmw[tmp], marker='.', s = 220, facecolors='none', edgecolors='blue', label=r'DECam LS dr8 in DESI')
if survey == 'north':
    tmp = (desitiles['in_desi'] ==  1) & (desitiles['dec'] > dec_resol_ns)
    SC         = ax.scatter(ramw[tmp],decmw[tmp], marker='.', s = 220, facecolors='none', edgecolors='blue', label=r'BASS/MzLS dr8 in DESI')


if survey == 'south':
    desitilesindes = get_isdes(desitiles['ra'],desitiles['dec'])
    tmp = (desitiles['in_desi'] ==  1) & (desitiles['dec'] < dec_resol_ns) & (desitilesindes)
    SC         = ax.scatter(ramw[tmp],decmw[tmp], marker='.', s = 220, facecolors='none', edgecolors='red', label=r'DECam DES in DESI')

gc = SkyCoord(l=np.linspace(0, 360, 50)*units.degree, b=np.full(50, 0)*units.degree, frame='galactic')
ramw,decmw = get_radec_mw(gc.icrs.ra.value, gc.icrs.dec.value,org)
SC         = ax.scatter(ramw, decmw, color='r', marker='.', lw=2)

lgnd = ax.legend()
[handle.set_sizes([120.0]) for handle in lgnd.legendHandles]

fig.savefig('%s/skyplot_%s.png' %(pathdir, survey), bbox_inches = 'tight', pad_inches = 0)

area_all: 	 20332.475(whole) 	 14265.679(desifootprint)
bgsarea_south: 	 15173.940(whole) 	 9716.980(desifootprint)
bgsarea_north: 	 5158.535(whole) 	 4548.699(desifootprint)
bgsarea_des: 	 5060.633(whole) 	 1113.520(desifootprint)
NGC DECaLS in DESI footprint: 	 5322.602
SGC DECaLS in DESI footprint: 	 4394.378

get area by number of passes + skymap of regions

In [7]:
# get the area covered by the number of passes from 1 to 3 and per band and for the joint bands.
Apasses = {}
for i in [0, 1, 2]:
    for j in ['G', 'R', 'Z', 'all']:
        
        if j == 'all': mask = (ran['NOBS_G'] > i) & (ran['NOBS_R'] > i) & (ran['NOBS_Z'] > i)
        else: mask = (ran['NOBS_'+j] > i)
            
        hpdicttest = get_dict(pixmapfile=dr8pix, hppix_ran=hppix_ran, maskrand=mask, Nranfiles=Nranfiles, 
                              ranindesi=ranindesi, desifootprint=True, target_outputs=False)
        
        for k in ['south', 'north', 'decals', 'des']:
            Apasses[j+str(i+1)+'_'+k] = hpdicttest['bgsarea_'+k]
        Apasses[j+str(i+1)+'_'+'all'] = hpdicttest['area_all']
            
        print('%s >= %i: \t south: %2.3g, north: %2.3g, decals: %2.3g, des: %2.3g, all: %2.3g' 
              %(j, i+1, hpdicttest['bgsarea_south'], hpdicttest['bgsarea_north'], hpdicttest['bgsarea_decals'], 
                    hpdicttest['bgsarea_des'], hpdicttest['area_all']))

G >= 1: 	 south: 9.69e+03, north: 4.54e+03, decals: 8.57e+03, des: 1.11e+03, all: 1.42e+04
R >= 1: 	 south: 9.69e+03, north: 4.54e+03, decals: 8.57e+03, des: 1.11e+03, all: 1.42e+04
Z >= 1: 	 south: 9.69e+03, north: 4.54e+03, decals: 8.57e+03, des: 1.11e+03, all: 1.42e+04
all >= 1: 	 south: 9.67e+03, north: 4.52e+03, decals: 8.56e+03, des: 1.11e+03, all: 1.42e+04
G >= 2: 	 south: 9.45e+03, north: 4.49e+03, decals: 8.34e+03, des: 1.11e+03, all: 1.39e+04
R >= 2: 	 south: 9.42e+03, north: 4.49e+03, decals: 8.31e+03, des: 1.11e+03, all: 1.39e+04
Z >= 2: 	 south: 9.49e+03, north: 4.49e+03, decals: 8.37e+03, des: 1.11e+03, all: 1.4e+04
all >= 2: 	 south: 9.26e+03, north: 4.41e+03, decals: 8.15e+03, des: 1.11e+03, all: 1.37e+04
G >= 3: 	 south: 7.77e+03, north: 3.93e+03, decals: 6.66e+03, des: 1.11e+03, all: 1.17e+04
R >= 3: 	 south: 7.57e+03, north: 3.9e+03, decals: 6.46e+03, des: 1.11e+03, all: 1.15e+04
Z >= 3: 	 south: 8.04e+03, north: 3.89e+03, decals: 6.93e+03, des: 1.11e+03, all: 1.19e+04
all >= 3: 	 south: 6.87e+03, north: 3.18e+03, decals: 5.76e+03, des: 1.11e+03, all: 1e+04
In [8]:
#
print('BAND \t >= 1 \t\t >= 2 \t\t >= 3')
print('------------------------------------------')
for band in ['G', 'R', 'Z', 'all']:
    lab1 = band+'1'+'_south'
    lab2 = band+'2'+'_south'
    lab3 = band+'3'+'_south'
    
    print('%s: \t %.0f \t %.0f \t %.0f' %(band, Apasses[lab1], Apasses[lab2], Apasses[lab3]))

BAND 	 >= 1 		 >= 2 		 >= 3
------------------------------------------
G: 	 9687 	 9454 	 7769
R: 	 9686 	 9422 	 7569
Z: 	 9686 	 9487 	 8036
all: 	 9669 	 9257 	 6870

DEV ZONE

In [68]:
#bgslist = ['BS', 'LG', 'GC', 'nobs', 'FMC2', 'CC', 'QC_FM', 'QC_FI', 'QC_FF', 'QC_IVAR']
bgsbut20 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['QC_FM','QC_FI','QC_FF','QC_IVAR', 'FMC2', 'CC', 'nobs', 'BS', 'LG', 'GC'], bgsmask=bgsmask(), rlimit=20)
bgsbut201 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['QC_FM','QC_FI','QC_FF','QC_IVAR', 'FMC2', 'CC', 'nobs', 'BS', 'LG', 'GC'], bgsmask=bgsmask(), rlimit=20.1)
bgsbut205 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['QC_FM','QC_FI','QC_FF','QC_IVAR', 'FMC2', 'CC', 'nobs', 'BS', 'LG', 'GC'], bgsmask=bgsmask(), rlimit=20.5)
bgsbut207 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['QC_FM','QC_FI','QC_FF','QC_IVAR', 'FMC2', 'CC', 'nobs', 'BS', 'LG', 'GC'], bgsmask=bgsmask(), rlimit=20.7)
#bgsbutCC = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['CC'], bgsmask=bgsmask())
#bgsbutFMC2 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['FMC2'], bgsmask=bgsmask())
#bgsbutBSLG = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=['BS', 'LG'], bgsmask=bgsmask())
In [43]:
bgs_any = ((cat['BGSBITS'] & 2**(20)) != 0)
In [69]:
namesels = {}
#namesels['bgs_any'] = bgs_any
namesels['bgsbut20'] = bgsbut20
namesels['bgsbut201'] = bgsbut201
namesels['bgsbut205'] = bgsbut205
namesels['bgsbut207'] = bgsbut207
In [70]:
#dic with default BGS selection and in DESI footprint
#maskrand=((rancuts['BS']) & (rancuts['GC']) & (rancuts['LG']) & (rancuts['nobs']))
hpdict_tmp = get_dict(cat=cat, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=hppix_cat, 
                      maskrand=None,
                          maskcat=None, 
                             Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None, 
                                desifootprint=True, namesels=namesels, galb=galb, log=True, survey='custom')

positions and desifotprint DONE...
systematics DONE...
randdens =  15000  ; len randoms =  145754705
bgsfracarea DONE...
regions DONE...
areas DONE...
computing for  north / bgsbut20
computing for  north / bgsbut201
computing for  north / bgsbut205
computing for  north / bgsbut207
target densities in north DONE...
computing for  south / bgsbut20
computing for  south / bgsbut201
computing for  south / bgsbut205
computing for  south / bgsbut207
target densities in south DONE...
meandens_bgsbut20_all = 1452 /deg2
meandens_bgsbut20_des = 1442 /deg2
meandens_bgsbut20_decals = 1453 /deg2
meandens_bgsbut20_north = nan /deg2
meandens_bgsbut20_south = 1452 /deg2
meandens_bgsbut20_svfields = 1466 /deg2
meandens_bgsbut20_svfields_n = nan /deg2
meandens_bgsbut20_svfields_s = 1466 /deg2
meandens_bgsbut201_all = 1608 /deg2
meandens_bgsbut201_des = 1596 /deg2
meandens_bgsbut201_decals = 1610 /deg2
meandens_bgsbut201_north = nan /deg2
meandens_bgsbut201_south = 1608 /deg2
meandens_bgsbut201_svfields = 1625 /deg2
meandens_bgsbut201_svfields_n = nan /deg2
meandens_bgsbut201_svfields_s = 1625 /deg2
meandens_bgsbut205_all = 2405 /deg2
meandens_bgsbut205_des = 2377 /deg2
meandens_bgsbut205_decals = 2408 /deg2
meandens_bgsbut205_north = nan /deg2
meandens_bgsbut205_south = 2405 /deg2
meandens_bgsbut205_svfields = 2432 /deg2
meandens_bgsbut205_svfields_n = nan /deg2
meandens_bgsbut205_svfields_s = 2432 /deg2
meandens_bgsbut207_all = 2948 /deg2
meandens_bgsbut207_des = 2901 /deg2
meandens_bgsbut207_decals = 2954 /deg2
meandens_bgsbut207_north = nan /deg2
meandens_bgsbut207_south = 2948 /deg2
meandens_bgsbut207_svfields = 2972 /deg2
meandens_bgsbut207_svfields_n = nan /deg2
meandens_bgsbut207_svfields_s = 2972 /deg2
dens_bgsbut20_all = 1451 /deg2
dens_bgsbut20_des = 1441 /deg2
dens_bgsbut20_decals = 1452 /deg2
dens_bgsbut20_north = nan /deg2
dens_bgsbut20_south = 1451 /deg2
dens_bgsbut20_svfields = 1464 /deg2
dens_bgsbut20_svfields_n = nan /deg2
dens_bgsbut20_svfields_s = 1464 /deg2
dens_bgsbut201_all = 1607 /deg2
dens_bgsbut201_des = 1594 /deg2
dens_bgsbut201_decals = 1609 /deg2
dens_bgsbut201_north = nan /deg2
dens_bgsbut201_south = 1607 /deg2
dens_bgsbut201_svfields = 1623 /deg2
dens_bgsbut201_svfields_n = nan /deg2
dens_bgsbut201_svfields_s = 1623 /deg2
dens_bgsbut205_all = 2403 /deg2
dens_bgsbut205_des = 2374 /deg2
dens_bgsbut205_decals = 2407 /deg2
dens_bgsbut205_north = nan /deg2
dens_bgsbut205_south = 2403 /deg2
dens_bgsbut205_svfields = 2428 /deg2
dens_bgsbut205_svfields_n = nan /deg2
dens_bgsbut205_svfields_s = 2428 /deg2
dens_bgsbut207_all = 2946 /deg2
dens_bgsbut207_des = 2898 /deg2
dens_bgsbut207_decals = 2952 /deg2
dens_bgsbut207_north = nan /deg2
dens_bgsbut207_south = 2946 /deg2
dens_bgsbut207_svfields = 2968 /deg2
dens_bgsbut207_svfields_n = nan /deg2
dens_bgsbut207_svfields_s = 2968 /deg2
In [71]:
fig    = plt.figure(figsize=(18,5))
gs     = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)
isdesi     = (hpdict_tmp['isdesi']) & (hpdict_tmp['bgsfracarea']>0)

plot_sysdens(hpdicttmp=hpdict_tmp, namesels=['bgsbut20', 'bgsbut201', 'bgsbut205', 'bgsbut207'], regs=[reg], 
                 syst='log10_stardens', mainreg=isdesi, xlim=None, n=0, nx=20, clip=True, 
                         denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test', label=True)

#plot_sysdens(hpdicttmp=hpdict_tmp, namesels=[ 'bgs_19.0_19.5', 'bgs_19.5_20.0'], regs=[reg], syst='log10_stardens', mainreg=isdesi, xlim=None, n=1, nx=nx, clip=True, 
#                         denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test', label=True)
In [73]:
#bgslist = ['BS', 'LG', 'GC', 'nobs', 'FMC2', 'CC', 'QC_FM', 'QC_FI', 'QC_FF', 'QC_IVAR']
bgsbut20 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=None, bgsmask=bgsmask(), rlimit=20)
bgsbut201 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=None, bgsmask=bgsmask(), rlimit=20.1)
bgsbut205 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=None, bgsmask=bgsmask(), rlimit=20.5)
bgsbut207 = bgsbut(bgsbits=cat['BGSBITS'], rmag=cat['RMAG'], pop=None, bgsmask=bgsmask(), rlimit=20.7)
In [74]:
namesels = {}
namesels['bgsbut20'] = bgsbut20
namesels['bgsbut201'] = bgsbut201
namesels['bgsbut205'] = bgsbut205
namesels['bgsbut207'] = bgsbut207
In [75]:
#dic with default BGS selection and in DESI footprint
#maskrand=((rancuts['BS']) & (rancuts['GC']) & (rancuts['LG']) & (rancuts['nobs']))
hpdict_tmp = get_dict(cat=cat, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=hppix_cat, 
                      maskrand=((rancuts['BS']) & (rancuts['GC']) & (rancuts['LG']) & (rancuts['nobs'])),
                          maskcat=None, 
                             Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None, 
                                desifootprint=True, namesels=namesels, galb=galb, log=True, survey='custom')

positions and desifotprint DONE...
systematics DONE...
randdens =  15000  ; len randoms =  145754705
bgsfracarea DONE...
regions DONE...
areas DONE...
computing for  north / bgsbut20
computing for  north / bgsbut201
computing for  north / bgsbut205
computing for  north / bgsbut207
target densities in north DONE...
computing for  south / bgsbut20
computing for  south / bgsbut201
computing for  south / bgsbut205
computing for  south / bgsbut207
target densities in south DONE...
meandens_bgsbut20_all = 1426 /deg2
meandens_bgsbut20_des = 1430 /deg2
meandens_bgsbut20_decals = 1425 /deg2
meandens_bgsbut20_north = nan /deg2
meandens_bgsbut20_south = 1426 /deg2
meandens_bgsbut20_svfields = 1437 /deg2
meandens_bgsbut20_svfields_n = nan /deg2
meandens_bgsbut20_svfields_s = 1437 /deg2
meandens_bgsbut201_all = 1578 /deg2
meandens_bgsbut201_des = 1583 /deg2
meandens_bgsbut201_decals = 1578 /deg2
meandens_bgsbut201_north = nan /deg2
meandens_bgsbut201_south = 1578 /deg2
meandens_bgsbut201_svfields = 1593 /deg2
meandens_bgsbut201_svfields_n = nan /deg2
meandens_bgsbut201_svfields_s = 1593 /deg2
meandens_bgsbut205_all = 2357 /deg2
meandens_bgsbut205_des = 2356 /deg2
meandens_bgsbut205_decals = 2357 /deg2
meandens_bgsbut205_north = nan /deg2
meandens_bgsbut205_south = 2357 /deg2
meandens_bgsbut205_svfields = 2384 /deg2
meandens_bgsbut205_svfields_n = nan /deg2
meandens_bgsbut205_svfields_s = 2384 /deg2
meandens_bgsbut207_all = 2890 /deg2
meandens_bgsbut207_des = 2876 /deg2
meandens_bgsbut207_decals = 2892 /deg2
meandens_bgsbut207_north = nan /deg2
meandens_bgsbut207_south = 2890 /deg2
meandens_bgsbut207_svfields = 2915 /deg2
meandens_bgsbut207_svfields_n = nan /deg2
meandens_bgsbut207_svfields_s = 2915 /deg2
dens_bgsbut20_all = 1423 /deg2
dens_bgsbut20_des = 1428 /deg2
dens_bgsbut20_decals = 1423 /deg2
dens_bgsbut20_north = nan /deg2
dens_bgsbut20_south = 1423 /deg2
dens_bgsbut20_svfields = 1435 /deg2
dens_bgsbut20_svfields_n = nan /deg2
dens_bgsbut20_svfields_s = 1435 /deg2
dens_bgsbut201_all = 1576 /deg2
dens_bgsbut201_des = 1581 /deg2
dens_bgsbut201_decals = 1575 /deg2
dens_bgsbut201_north = nan /deg2
dens_bgsbut201_south = 1576 /deg2
dens_bgsbut201_svfields = 1591 /deg2
dens_bgsbut201_svfields_n = nan /deg2
dens_bgsbut201_svfields_s = 1591 /deg2
dens_bgsbut205_all = 2354 /deg2
dens_bgsbut205_des = 2353 /deg2
dens_bgsbut205_decals = 2354 /deg2
dens_bgsbut205_north = nan /deg2
dens_bgsbut205_south = 2354 /deg2
dens_bgsbut205_svfields = 2381 /deg2
dens_bgsbut205_svfields_n = nan /deg2
dens_bgsbut205_svfields_s = 2381 /deg2
dens_bgsbut207_all = 2885 /deg2
dens_bgsbut207_des = 2872 /deg2
dens_bgsbut207_decals = 2887 /deg2
dens_bgsbut207_north = nan /deg2
dens_bgsbut207_south = 2885 /deg2
dens_bgsbut207_svfields = 2912 /deg2
dens_bgsbut207_svfields_n = nan /deg2
dens_bgsbut207_svfields_s = 2912 /deg2
In [76]:
fig    = plt.figure(figsize=(18,5))
gs     = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)
isdesi     = (hpdict_tmp['isdesi']) & (hpdict_tmp['bgsfracarea']>0)

plot_sysdens(hpdicttmp=hpdict_tmp, namesels=['bgsbut20', 'bgsbut201', 'bgsbut205', 'bgsbut207'], regs=[reg], 
                 syst='log10_stardens', mainreg=isdesi, xlim=None, n=0, nx=20, clip=True, 
                         denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test2', label=True)

#plot_sysdens(hpdicttmp=hpdict_tmp, namesels=[ 'bgs_19.0_19.5', 'bgs_19.5_20.0'], regs=[reg], syst='log10_stardens', mainreg=isdesi, xlim=None, n=1, nx=nx, clip=True, 
#                         denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test', label=True)
In [80]:
#
ws = 1./((-3.96*10**(-5))*hpdict['stardens'] + 1.03)

fig    = plt.figure(figsize=(18,5))
gs     = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)

plot_sysdens(hpdicttmp=hpdict, namesels=['bgs_any', 'bright', 'faint'], regs=[reg], syst='log10_stardens', mainreg=isdesi, xlim=None, n=0, nx=nx, clip=True, 
                         denslims=False, ylab=True, weights=True, fig=fig, gs=gs, title='test', label=True, ws=ws, onlyweights=True)
In [88]:
hpdict['meandens_bgs_any_south']
Out[88]:
1425.646415390457
In [89]:
hpdict_tmp['meandens_bgs_any_south']
Out[89]:
1432.660001909052
In [90]:
#
#ws = 1./((-3.96*10**(-5))*hpdict['stardens'] + 1.03)

fig    = plt.figure(figsize=(18,5))
gs     = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)

plot_sysdens(hpdicttmp=hpdict_tmp, namesels=['bgs_any', 'bright', 'faint'], regs=[reg], syst='log10_stardens', mainreg=isdesi, xlim=None, n=0, nx=nx, clip=True, 
                         denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test', label=True, ws=None, onlyweights=False)
In [53]:
bgs_any = ((cat['BGSBITS'] & 2**(20)) != 0)
namesels = {}
mags = [15, 18, 18.5, 19.0, 19.5, 20.0]
for num,i in enumerate(mags[:-1]):
    namesels['bgs_%s_%s' %(str(mags[num]), str(mags[num+1]))] = (bgs_any) & (cat['RMAG'] > mags[num]) & (cat['RMAG'] < mags[num+1])
In [55]:
#dic with default BGS selection and in DESI footprint
hpdict_tmp = get_dict(cat=cat, pixmapfile=dr8pix, hppix_ran=hppix_ran, hppix_cat=hppix_cat, 
                      maskrand=((rancuts['BS']) & (rancuts['GC']) & (rancuts['LG']) & (rancuts['nobs'])),
                          maskcat=None, 
                             Nranfiles=Nranfiles, ranindesi=ranindesi, catindesi=None, 
                                desifootprint=True, namesels=namesels, galb=galb, log=True, survey='custom')

positions and desifotprint DONE...
systematics DONE...
randdens =  15000  ; len randoms =  145754705
bgsfracarea DONE...
regions DONE...
areas DONE...
computing for  north / bgs_15_18
computing for  north / bgs_18_18.5
computing for  north / bgs_18.5_18.5
computing for  north / bgs_18.5_19.0
computing for  north / bgs_19.0_19.5
computing for  north / bgs_19.5_20.0
target densities in north DONE...
computing for  south / bgs_15_18
computing for  south / bgs_18_18.5
computing for  south / bgs_18.5_18.5
computing for  south / bgs_18.5_19.0
computing for  south / bgs_19.0_19.5
computing for  south / bgs_19.5_20.0
target densities in south DONE...
meandens_bgs_15_18_all = 154 /deg2
meandens_bgs_15_18_des = 153 /deg2
meandens_bgs_15_18_decals = 154 /deg2
meandens_bgs_15_18_north = nan /deg2
meandens_bgs_15_18_south = 154 /deg2
meandens_bgs_15_18_svfields = 151 /deg2
meandens_bgs_15_18_svfields_n = nan /deg2
meandens_bgs_15_18_svfields_s = 151 /deg2
meandens_bgs_18_18.5_all = 124 /deg2
meandens_bgs_18_18.5_des = 124 /deg2
meandens_bgs_18_18.5_decals = 124 /deg2
meandens_bgs_18_18.5_north = nan /deg2
meandens_bgs_18_18.5_south = 124 /deg2
meandens_bgs_18_18.5_svfields = 123 /deg2
meandens_bgs_18_18.5_svfields_n = nan /deg2
meandens_bgs_18_18.5_svfields_s = 123 /deg2
meandens_bgs_18.5_18.5_all = 0 /deg2
meandens_bgs_18.5_18.5_des = 0 /deg2
meandens_bgs_18.5_18.5_decals = 0 /deg2
meandens_bgs_18.5_18.5_north = nan /deg2
meandens_bgs_18.5_18.5_south = 0 /deg2
meandens_bgs_18.5_18.5_svfields = 0 /deg2
meandens_bgs_18.5_18.5_svfields_n = nan /deg2
meandens_bgs_18.5_18.5_svfields_s = 0 /deg2
meandens_bgs_18.5_19.0_all = 211 /deg2
meandens_bgs_18.5_19.0_des = 213 /deg2
meandens_bgs_18.5_19.0_decals = 211 /deg2
meandens_bgs_18.5_19.0_north = nan /deg2
meandens_bgs_18.5_19.0_south = 211 /deg2
meandens_bgs_18.5_19.0_svfields = 212 /deg2
meandens_bgs_18.5_19.0_svfields_n = nan /deg2
meandens_bgs_18.5_19.0_svfields_s = 212 /deg2
meandens_bgs_19.0_19.5_all = 354 /deg2
meandens_bgs_19.0_19.5_des = 356 /deg2
meandens_bgs_19.0_19.5_decals = 353 /deg2
meandens_bgs_19.0_19.5_north = nan /deg2
meandens_bgs_19.0_19.5_south = 354 /deg2
meandens_bgs_19.0_19.5_svfields = 358 /deg2
meandens_bgs_19.0_19.5_svfields_n = nan /deg2
meandens_bgs_19.0_19.5_svfields_s = 358 /deg2
meandens_bgs_19.5_20.0_all = 579 /deg2
meandens_bgs_19.5_20.0_des = 580 /deg2
meandens_bgs_19.5_20.0_decals = 578 /deg2
meandens_bgs_19.5_20.0_north = nan /deg2
meandens_bgs_19.5_20.0_south = 579 /deg2
meandens_bgs_19.5_20.0_svfields = 589 /deg2
meandens_bgs_19.5_20.0_svfields_n = nan /deg2
meandens_bgs_19.5_20.0_svfields_s = 589 /deg2
dens_bgs_15_18_all = 154 /deg2
dens_bgs_15_18_des = 153 /deg2
dens_bgs_15_18_decals = 154 /deg2
dens_bgs_15_18_north = nan /deg2
dens_bgs_15_18_south = 154 /deg2
dens_bgs_15_18_svfields = 151 /deg2
dens_bgs_15_18_svfields_n = nan /deg2
dens_bgs_15_18_svfields_s = 151 /deg2
dens_bgs_18_18.5_all = 124 /deg2
dens_bgs_18_18.5_des = 124 /deg2
dens_bgs_18_18.5_decals = 124 /deg2
dens_bgs_18_18.5_north = nan /deg2
dens_bgs_18_18.5_south = 124 /deg2
dens_bgs_18_18.5_svfields = 122 /deg2
dens_bgs_18_18.5_svfields_n = nan /deg2
dens_bgs_18_18.5_svfields_s = 122 /deg2
dens_bgs_18.5_18.5_all = 0 /deg2
dens_bgs_18.5_18.5_des = 0 /deg2
dens_bgs_18.5_18.5_decals = 0 /deg2
dens_bgs_18.5_18.5_north = nan /deg2
dens_bgs_18.5_18.5_south = 0 /deg2
dens_bgs_18.5_18.5_svfields = 0 /deg2
dens_bgs_18.5_18.5_svfields_n = nan /deg2
dens_bgs_18.5_18.5_svfields_s = 0 /deg2
dens_bgs_18.5_19.0_all = 211 /deg2
dens_bgs_18.5_19.0_des = 213 /deg2
dens_bgs_18.5_19.0_decals = 211 /deg2
dens_bgs_18.5_19.0_north = nan /deg2
dens_bgs_18.5_19.0_south = 211 /deg2
dens_bgs_18.5_19.0_svfields = 212 /deg2
dens_bgs_18.5_19.0_svfields_n = nan /deg2
dens_bgs_18.5_19.0_svfields_s = 212 /deg2
dens_bgs_19.0_19.5_all = 353 /deg2
dens_bgs_19.0_19.5_des = 356 /deg2
dens_bgs_19.0_19.5_decals = 353 /deg2
dens_bgs_19.0_19.5_north = nan /deg2
dens_bgs_19.0_19.5_south = 353 /deg2
dens_bgs_19.0_19.5_svfields = 357 /deg2
dens_bgs_19.0_19.5_svfields_n = nan /deg2
dens_bgs_19.0_19.5_svfields_s = 357 /deg2
dens_bgs_19.5_20.0_all = 578 /deg2
dens_bgs_19.5_20.0_des = 579 /deg2
dens_bgs_19.5_20.0_decals = 578 /deg2
dens_bgs_19.5_20.0_north = nan /deg2
dens_bgs_19.5_20.0_south = 578 /deg2
dens_bgs_19.5_20.0_svfields = 588 /deg2
dens_bgs_19.5_20.0_svfields_n = nan /deg2
dens_bgs_19.5_20.0_svfields_s = 588 /deg2
In [60]:
fig    = plt.figure(figsize=(18,5))
gs     = gridspec.GridSpec(1,2,hspace=0.30,wspace=0.25)

plot_sysdens(hpdicttmp=hpdict_tmp, namesels=['bgs_15_18', 'bgs_18_18.5', 'bgs_18.5_19.0'], regs=[reg], syst='log10_stardens', mainreg=isdesi, xlim=None, n=0, nx=nx, clip=True, 
                         denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test', label=True)

plot_sysdens(hpdicttmp=hpdict_tmp, namesels=[ 'bgs_19.0_19.5', 'bgs_19.5_20.0'], regs=[reg], syst='log10_stardens', mainreg=isdesi, xlim=None, n=1, nx=nx, clip=True, 
                         denslims=False, ylab=True, weights=False, fig=fig, gs=gs, title='test', label=True)
In [ ]:
bgs_bright = ((cat['BGSBITS'] & 2**(21)) != 0)
bgs_faint = ((cat['BGSBITS'] & 2**(22)) != 0)
bgs_any = ((cat['BGSBITS'] & 2**(20)) != 0)